
params <- matrix(nrow=1080,ncol=6)
dpl <- c(0,0.25,0.5,0.75,1)
scl <- c(0.001,0.01,0.1,0.5,1,2)
sml <- c(10^-5,10^-6,10^-7,10^-8)
nml <- c(10^-4,10^-5,10^-7)
sw <- 0

for (divprob in 1:5){
  for (selcoeff in 1:6){
    for (selmu in 1:4){
      for (neutmu in 1:3){
        for (replicates in 1:3){
          sw <- sw + 1
          params[sw,1] <- sw
          params[sw,2] <- dpl[divprob]
          params[sw,3] <- scl[selcoeff]
          params[sw,4] <- sml[selmu]
          params[sw,5] <- nml[neutmu]
          params[sw,6] <- replicates
        }
      }
    }	
  }
}

# Choose which sweeps to plot in whittaker
sweeps <- params[params[,3]==2 & params[,4]==10^-7 & params[,5]==10^-5,1]

pdf(paste("whittaker-plot.pdf",sep=""),width=6,height=6,pointsize=12)
for(ii in 1:length(sweeps)){
  i <- sweeps[ii]
  sim.title <- paste("sim=",params[i,1]," r=",params[i,2]," sm=",params[i,4]," nm=",params[i,5]," s=",params[i,3]," run=",params[i,6],sep="")
  grid.log <- read.csv(paste("../sim",i,"/sim",i,".grid.log",sep=""),header=FALSE)
  clones.to.patterns <- read.csv(paste("../sim",i,"/sim",i,".clones.to.patterns",sep=""),header=FALSE)

  # take end of simulaiton
  grid.log <- grid.log[grid.log[,1]==7300,]

  clones <- grid.log[grid.log[,4]!=-1,4]
  clone.counts <- as.numeric(sort(table(clones),decreasing=TRUE))
  total.clones <- sum(clone.counts)
  alls <- length(clone.counts) 
  
  # Plot the top 100 species
  plot(c(1:100),log10(clone.counts[1:100]/total.clones),ylim=log10(c(0.0001,1)),xlim=c(1,100),axes=F,xlab="Species rank",ylab="Relative abundance",pch=".")
  title(main=sim.title)
  axis(2,at=log10(c(0.0001,0.001,0.01,0.1,1)),label=c(0.0001,0.001,0.01,0.1,1))
  axis(1,at=c(1,25,50,75,100),label=c(1,25,50,75,100))

  # Plot all species  
plot(c(1:alls),log10(clone.counts[1:alls]/total.clones),ylim=log10(c(0.00001,1)),xlim=c(1,alls),axes=F,xlab="Species rank",ylab="Relative abundance",pch=".")
  title(main=sim.title)
  axis(2,at=log10(c(0.00001,0.0001,0.001,0.01,0.1,1)),label=c(0.00001,0.0001,0.001,0.01,0.1,1))
  axis(1,at=c(1,round(alls/2),alls),label=c(1,round(alls/2),alls))
  
  # Histogram of top 100 species
  hist(clone.counts[1:100],breaks=50,xlab="Number of individuals",ylab="Number of species",freq=TRUE,main=NULL)
  title(main=sim.title)
    
  # Histogram of all species
  hist(clone.counts[1:alls],breaks=50,xlab="Number of individuals",ylab="Number of species",freq=TRUE,main=NULL)
  title(main=sim.title)
  print(paste("sim",i,"done."))
}

dev.off()
